Memory Reactivation & Update

Memory Reactivation & Update

Deep Memory Lab, Cognitive Science Research Group, Korea Brain Research Institute

1 Outline

2 Results

set.seed(12345) # for reproducibility
options(knitr.kable.NA = '')

# install // load packages 
# Some packages need to be loaded. 
# We use `pacman` as a package manager, which takes care of the other packages. 
if (!require("distill", quietly = TRUE)) install.packages("distill")
if (!require("devtools", quietly = TRUE)) install.packages("devtools")
if (!require("papaja", quietly = TRUE)) devtools::install_github("crsh/papaja")
if (!require("patchwork", quietly = TRUE)) devtools::install_github("thomasp85/patchwork")
if (!require("klippy", quietly = TRUE)) devtools::install_github("RLesur/klippy")
if (!require("pacman", quietly = TRUE)) install.packages("pacman")
if (!require("Rmisc", quietly = TRUE)) install.packages("Rmisc")
if (!require("rstatix", quietly = TRUE)) install.packages("rstatix")
if (!require("effsize", quietly = TRUE)) install.packages("effsize")
if (!require("lsr", quietly = TRUE)) install.packages("lsr")
if (!require("effectsize", quietly = TRUE)) install.packages("effectsize")
if (!require("ggbeeswarm", quietly = TRUE)) install.packages("ggbeeswarm") # Never load it directly.
pacman::p_load(tidyverse, papaja, knitr, dplyr, car, psych, afex, lme4, lmerTest, 
               emmeans, ggplot2, ggpubr, lattice, latticeExtra, parallel, effects, psycho, apa)
library("patchwork"); library("klippy")
klippy::klippy()

2.1 Phase 1 - Categorization

t1 <- read.csv("data/reactMEM_p1_all.csv", header = T)
glimpse(t1, width = 70)
## Rows: 32,640
## Columns: 15
## $ SN      <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ Trial   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
## $ Block   <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ bTrial  <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
## $ cCong   <int> 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 2, 2, 2, 2, 1, 1, 2, 2…
## $ cnRep   <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2…
## $ cnPair  <int> 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2…
## $ cIdt    <int> 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 4, 4, 5, 5, 3, 3…
## $ cCat    <int> 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 1, 2, 1, 2, 1…
## $ Im_nExm <int> 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1…
## $ Im_Idx  <int> 199, 193, 47, 77, 157, 127, 209, 111, 181, 72, 61, 1…
## $ Im_Name <chr> "199_3_out.jpg", "193_1_out.jpg", "047_3_in.jpg", "0…
## $ Resp    <int> 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 1, 2, 1, 2, 1…
## $ RT      <dbl> 1.1000, 1.0478, 0.7386, 0.7283, 0.8903, 0.9315, 0.70…
## $ Corr    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
length(unique(t1$SN))
## [1] 32

# change class of main factors: double to factor
t1$SN = factor(t1$SN)
t1$Block = factor(t1$Block)
t1$cCong = factor(t1$cCong, levels=c(1,2), labels=c("v","nv")) 
t1$cnRep = factor(t1$cnRep, levels=c(1,2,3,5,4), labels=c(1,2,3,4,0))
t1$cnPair = factor(t1$cnPair)
t1$cIdt = factor(t1$cIdt)
t1$cCat = factor(t1$cCat)
t1$Im_Name = factor(t1$Im_Name)
t1$Im_Idx = factor(t1$Im_Idx)
t1$RT <- t1$RT*1000;
t1$Corr <- as.numeric(t1$Corr==1)

# check number of trials for each condition/SN
table(t1$SN)
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
## 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 
##   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32 
## 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020 1020
table(t1$Block, t1$SN) 
##    
##       1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19
##   1 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170
##   2 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170
##   3 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170
##   4 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170
##   5 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170
##   6 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170 170
##    
##      20  21  22  23  24  25  26  27  28  29  30  31  32
##   1 170 170 170 170 170 170 170 170 170 170 170 170 170
##   2 170 170 170 170 170 170 170 170 170 170 170 170 170
##   3 170 170 170 170 170 170 170 170 170 170 170 170 170
##   4 170 170 170 170 170 170 170 170 170 170 170 170 170
##   5 170 170 170 170 170 170 170 170 170 170 170 170 170
##   6 170 170 170 170 170 170 170 170 170 170 170 170 170
table(t1$cCong, t1$SN) 
##     
##        1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
##   v  540 540 540 540 540 540 540 540 540 540 540 540 540 540 540 540 540 540
##   nv 480 480 480 480 480 480 480 480 480 480 480 480 480 480 480 480 480 480
##     
##       19  20  21  22  23  24  25  26  27  28  29  30  31  32
##   v  540 540 540 540 540 540 540 540 540 540 540 540 540 540
##   nv 480 480 480 480 480 480 480 480 480 480 480 480 480 480
headTail(t1)
##         SN Trial Block bTrial cCong cnRep cnPair cIdt cCat Im_nExm Im_Idx
## 1        1     1     1      1     v     1      1    1    2       3    199
## 2        1     2     1      2     v     1      2    1    2       1    193
## 3        1     3     1      3     v     1      1    2    1       3     47
## 4        1     4     1      4     v     1      2    2    1       1     77
## ...   <NA>   ...  <NA>    ...  <NA>  <NA>   <NA> <NA> <NA>     ...   <NA>
## 32637   32  1017     6    167    nv     0      2  118    1       3     54
## 32638   32  1018     6    168     v     4      2  119    2       1    144
## 32639   32  1019     6    169    nv     4      2  120    1       1     13
## 32640   32  1020     6    170    nv     4      2  118    1       1     54
##             Im_Name Resp     RT Corr
## 1     199_3_out.jpg    2   1100    1
## 2     193_1_out.jpg    2 1047.8    1
## 3      047_3_in.jpg    1  738.6    1
## 4      077_1_in.jpg    1  728.3    1
## ...            <NA>  ...    ...  ...
## 32637  054_3_in.jpg    1 1392.9    1
## 32638 144_1_out.jpg    2  765.9    1
## 32639  013_1_in.jpg    1 1465.3    1
## 32640  054_1_in.jpg    1 1380.5    1
glimpse(t1, width = 70)
## Rows: 32,640
## Columns: 15
## $ SN      <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ Trial   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
## $ Block   <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ bTrial  <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
## $ cCong   <fct> v, v, v, v, nv, nv, nv, nv, v, v, nv, nv, nv, nv, v,…
## $ cnRep   <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2…
## $ cnPair  <fct> 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2…
## $ cIdt    <fct> 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 4, 4, 5, 5, 3, 3…
## $ cCat    <fct> 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 1, 2, 1, 2, 1…
## $ Im_nExm <int> 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1…
## $ Im_Idx  <fct> 199, 193, 47, 77, 157, 127, 209, 111, 181, 72, 61, 1…
## $ Im_Name <fct> 199_3_out.jpg, 193_1_out.jpg, 047_3_in.jpg, 077_1_in…
## $ Resp    <int> 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 1, 2, 1, 2, 1…
## $ RT      <dbl> 1100.0, 1047.8, 738.6, 728.3, 890.3, 931.5, 701.8, 8…
## $ Corr    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…

2.1.1 RT analysis: General

## 1.2.1. preprocessing
# filter null/incorrect data
ct1 <- t1 %>% filter(Corr == 1) # filter(Corr == 1) # remove incorrect trial

# check accuracy
100-100*(nrow(ct1)/nrow(t1))
## [1] 4.666054
# 5.76% incorrect trials were not analyzed

# trimming 3sd outlier trials
# tt1 <- ct1
tt1 <- ct1 %>% filter(RT > 200 & RT < 10000) %>%
  group_by(SN) %>% # grouping by participants
  nest() %>%
  mutate(lbound = map(data, ~mean(.$RT)-3*sd(.$RT)),
         ubound = map(data, ~mean(.$RT)+3*sd(.$RT))) %>% # make new data (3sd cut)
  unnest(c(lbound, ubound))%>%
  unnest(data) %>%
  mutate(Outlier = (RT < lbound)|(RT > ubound)) %>% # set outlier
  filter(Outlier == FALSE) %>% # filtering outlier
  ungroup() %>%
  dplyr::select(SN, Block, cCong, cnRep, cnPair, cIdt, cCat, Im_Name, RT, Corr) # select variables to analyze

# outlier trial ratio
100-100*(nrow(tt1)/nrow(ct1))
## [1] 1.812514
## [1] 1.89857

# RTs were trimmed
# when 1) faster than 200ms 2) slower than 10 secs 3) 3SD away from the participants mean

# mean number of trials for each conditions
tt1 %>% group_by(SN, cCong, cnRep) %>% 
  filter(cnPair == 2) %>% 
  summarise(NumTrial = length(RT)) %>%
  ungroup() %>%
  group_by(cCong, cnRep) %>%
  summarise(Mean = mean(NumTrial), 
            Median = median(NumTrial), 
            Min = min(NumTrial), 
            Max = max(NumTrial)) %>% 
  ungroup %>%
  kable(digits=2)
## `summarise()` has grouped output by 'SN', 'cCong'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'cCong'. You can override using the `.groups` argument.

cCong	cnRep	Mean	Median	Min	Max
v	1	55.06	56.0	48	58
v	2	56.00	56.0	51	59
v	3	56.44	58.0	50	59
v	4	55.97	56.0	44	60
v	0	56.19	56.5	49	60
nv	1	55.38	56.0	42	59
nv	2	56.22	57.0	51	60
nv	3	57.19	58.0	48	60
nv	4	57.22	57.0	52	60
nv	0	56.59	57.0	49	60

# check Distribution

# before trimming
den1 <- ggplot(ct1, aes(x=RT)) + 
  geom_density() + 
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) 

# after trimming
den2 <- ggplot(tt1, aes(x=RT)) + 
  geom_density() + 
  theme_bw(base_size = 18) + 
  labs(x = "Trimmed RT") + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) 
den1 + den2

ggplot(data=tt1, aes(x = SN, y = RT)) + 
  ggtitle("mean RT") +
  xlab("Participants") + ylab("RT") + theme_bw() +
  theme(text=element_text(size=14)) +
  theme(legend.title = element_text(face = 1,size = 15)) +
  geom_boxplot()

# subject-level, long format
t1rtL <- tt1 %>% filter(cnPair == 2) %>% group_by(SN, cCong, cnRep) %>%
  summarise(RT = mean(RT)) %>%
  ungroup()
## `summarise()` has grouped output by 'SN', 'cCong'. You can override using the `.groups` argument.
#t1rtL %>% kable(digits=2)

# subject-level, wide format (SN/Btw/Block)
t1rtW <- t1rtL %>% spread(key=cnRep, value = RT)
#t1rtW %>% kable(digits=2)

# subject-level, wide format (SN/Btw/Block)
t1rtW1 <- t1rtL %>% # filter(cnRep == 4) %>% 
  spread(key=cCong, value = RT)
#t1rtW1 %>% kable(digits=2)

# summary table: grand mean (eyerep/locrep)
t1rtG <- t1rtL %>% group_by(cCong, cnRep) %>%
  summarise(RT.m = mean(RT), RT.sd = sd(RT)) %>%
  ungroup()
## `summarise()` has grouped output by 'cCong'. You can override using the `.groups` argument.
t1rtG$RT.se <- Rmisc::summarySEwithin(data = t1rtL, measurevar = "RT", 
                                      idvar = "SN", withinvars = c("cCong", "cnRep"))$se
t1rtG$RT.ci <- Rmisc::summarySEwithin(data = t1rtL, measurevar = "RT", 
                                      idvar = "SN", withinvars = c("cCong", "cnRep"))$ci
t1rtG <- t1rtG %>% 
  mutate(lower.ci = RT.m-RT.ci,
         upper.ci = RT.m+RT.ci)
t1rtG %>% dplyr::select(cCong, cnRep, RT.m) %>% 
  spread(key=cnRep, value=RT.m) %>% kable(digits=2)

cCong	1	2	3	4	0
v	758.18	704.12	667.56	678.37	733.66
nv	760.14	700.95	677.13	663.65	734.82

p1.L.rt.g <- t1rtG %>% filter(cnRep != 0)
p1.L.rt.g$Cong <- factor(p1.L.rt.g$cCong, labels=c("Violation", "Non-violation"))
p1.L.rt.g.plot1 <- ggplot(data=p1.L.rt.g , 
                              aes(x=cnRep, y=RT.m, ymin=RT.m-RT.ci, ymax=RT.m+RT.ci, color=Cong, shape=Cong)) +
  geom_point(size = 5, position = position_dodge(.3)) +
  geom_errorbar(width = .2, position = position_dodge(.3)) +
  geom_line(aes(group = Cong), position = position_dodge(.3), size=1) + 
  scale_color_manual(values = c("#ED7D31", "#70AD47")) + # c("red", "black")) +
  # scale_shape_manual(values = c("red", "black"), labels = c("Experimental", "Control")) +
  coord_cartesian(ylim = c(500,1000), clip = "on") +
  # scale_y_continuous(breaks = seq(0,80)) +
  scale_x_discrete(labels = c("1", "2", "3", "4")) +
  labs(x = "Repetition", y = "Response Time (ms)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        legend.position = c(0.65, 0.85),
        legend.title = element_blank(),
        # aspect.ratio = 0.6,
        legend.background = element_blank(),
        plot.margin = margin(1, 0.3, 0.3, 0.5, "cm"),
        legend.key = element_blank())

p1.L.rt.g <- t1rtG %>% filter(cnRep != 2, cnRep != 3, cnRep != 4)
p1.L.rt.g$Cong <- factor(p1.L.rt.g$cCong, labels=c("Violation", "Non-violation"))
p1.L.rt.g.plot2 <- ggplot(data=p1.L.rt.g , 
                          aes(x=cnRep, y=RT.m, ymin=RT.m-RT.ci, ymax=RT.m+RT.ci, color=Cong, shape=Cong)) +
  geom_point(size = 5, position = position_dodge(.3)) +
  geom_errorbar(width = .2, position = position_dodge(.3)) +
  geom_line(aes(group = Cong), position = position_dodge(.3), size=1) + 
  scale_color_manual(values = c("#ED7D31", "#70AD47")) + # c("red", "black")) +
  # scale_shape_manual(values = c("red", "black"), labels = c("Experimental", "Control")) +
  coord_cartesian(ylim = c(500,1000), clip = "on") +
  # scale_y_continuous(breaks = seq(0,80)) +
  scale_x_discrete(labels = c("1st B", "B-Prime")) +
  labs(x = "Type", y = "Response Time (ms)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        legend.position = c(0.65, 0.85),
        legend.title = element_blank(),
        # aspect.ratio = 0.6,
        legend.background = element_blank(),
        plot.margin = margin(1, 0.3, 0.3, 0.5, "cm"),
        legend.key = element_blank())
# p1.L.rt.g.plot1
# p1.L.rt.g.plot2

p1.rt.plot <- ggarrange(p1.L.rt.g.plot1, p1.L.rt.g.plot2, ncol = 2,
          labels=c("A) Categorization of B", "B) Categorization fo B/B`"))
p1.rt.plot

RM ANOVA 1 - 반복에 따른 priming 효과

t1rtL.1 <- t1rtL %>% filter(cnRep != 0)
t1.rt.aov1 <- aov_ez(id="SN", dv = "RT", data = t1rtL.1, within = c("cCong", "cnRep"))
# test_sphericity(t2.acc.aov1)
# summary(t1.rt.aov1)
# anova(t1.rt.aov1, es = "pes")
nice(t1.rt.aov1, es="pes") %>% kable(digits=3)

Effect	df	MSE	F	pes	p.value
cCong	1, 31	869.19	0.19	.006	.670
cnRep	2.26, 70.16	1856.46	77.77 ***	.715	<.001
cCong:cnRep	2.32, 72.02	1261.63	1.70	.052	.185

RM ANOVA 1 - Post-hoc Test

### RM ANOVA 1 - post-hoc
t1.rt.m1 <- emmeans(t1.rt.aov1, pairwise ~ cCong | cnRep) # adjust="bon" , "tukey"
t1.rt.m1$contrasts %>% kable(digits=3)

contrast	cnRep	estimate	SE	df	t.ratio	p.value
v - nv	X1	-1.966	6.290	31	-0.313	0.757
v - nv	X2	3.165	7.914	31	0.400	0.692
v - nv	X3	-9.568	6.312	31	-1.516	0.140
v - nv	X4	14.717	9.772	31	1.506	0.142

plot(t1.rt.m1, horizontal = T, comparisons = T)

t1.rt.m2 <- emmeans(t1.rt.aov1, pairwise ~ cnRep | cCong) # adjust="bon" , "tukey"
t1.rt.m2$contrasts %>% kable(digits=3)

contrast	cCong	estimate	SE	df	t.ratio	p.value
X1 - X2	v	54.056	8.983	31	6.017	0.000
X1 - X3	v	90.618	9.704	31	9.338	0.000
X1 - X4	v	79.804	9.227	31	8.649	0.000
X2 - X3	v	36.563	8.987	31	4.068	0.002
X2 - X4	v	25.749	8.285	31	3.108	0.020
X3 - X4	v	-10.814	7.597	31	-1.423	0.495
X1 - X2	nv	59.187	7.354	31	8.048	0.000
X1 - X3	nv	83.016	7.993	31	10.386	0.000
X1 - X4	nv	96.487	10.198	31	9.462	0.000
X2 - X3	nv	23.829	7.141	31	3.337	0.011
X2 - X4	nv	37.300	9.800	31	3.806	0.003
X3 - X4	nv	13.471	7.448	31	1.809	0.289

plot(t1.rt.m2, horizontal = T, comparisons = T)

RM ANOVA 2 - violation 여부에 따른 priming 효과 감소 (differentiation effect)

t1rtL.2 <- t1rtL %>% filter(cnRep !=1, cnRep != 2, cnRep !=0)
t1.rt.aov1 <- aov_ez(id="SN", dv = "RT", data = t1rtL.2, within = c("cCong", "cnRep"))
# test_sphericity(t2.acc.aov1)
#summary(t1.rt.aov1)
#anova(t1.rt.aov1, es = "pes")
nice(t1.rt.aov1, es="pes") %>% kable(digits=3)

Effect	df	MSE	F	pes	p.value
cCong	1, 31	969.52	0.22	.007	.643
cnRep	1, 31	615.18	0.09	.003	.764
cCong:cnRep	1, 31	1195.72	3.95 +	.113	.056

RM ANOVA 2 - Post-hoc Test

### RM ANOVA 2 - post-hoc
t1.rt.m1 <- emmeans(t1.rt.aov1, pairwise ~ cnRep | cCong) # adjust="bon" , "tukey"
t1.rt.m1$contrasts %>% kable(digits=3)

contrast	cCong	estimate	SE	df	t.ratio	p.value
X3 - X4	v	-10.814	7.597	31	-1.423	0.165
X3 - X4	nv	13.471	7.448	31	1.809	0.080

plot(t1.rt.m1, horizontal = FALSE, comparisons = T)

t1.rt.m1 <- emmeans(t1.rt.aov1, pairwise ~ cCong | cnRep) # adjust="bon" , "tukey"
t1.rt.m1$contrasts %>% kable(digits=3)

contrast	cnRep	estimate	SE	df	t.ratio	p.value
v - nv	X3	-9.568	6.312	31	-1.516	0.140
v - nv	X4	14.717	9.772	31	1.506	0.142

plot(t1.rt.m1, horizontal = FALSE, comparisons = T)

RM ANOVA 3 - general category priming 효과

t1rtL.3 <- t1rtL %>% filter(cnRep !=2, cnRep !=3, cnRep !=4)
t1.rt.aov1 <- aov_ez(id="SN", dv = "RT", data = t1rtL.3, within = c("cCong", "cnRep"))
# test_sphericity(t2.acc.aov1)
# summary(t1.rt.aov1)
# anova(t1.rt.aov1, es = "pes")
nice(t1.rt.aov1, es="pes") %>% kable(digits=3)

Effect	df	MSE	F	pes	p.value
cCong	1, 31	883.47	0.09	.003	.768
cnRep	1, 31	1031.18	19.27 ***	.383	<.001
cCong:cnRep	1, 31	677.30	0.01	<.001	.931

RM ANOVA 3 - Post-hoc Test

t1.rt.m1 <- emmeans(t1.rt.aov1, pairwise ~ cnRep | cCong) # adjust="bon" , "tukey"
t1.rt.m1$contrasts %>% kable(digits=3)

contrast	cCong	estimate	SE	df	t.ratio	p.value
X1 - X0	v	24.518	8.743	31	2.804	0.009
X1 - X0	nv	25.320	5.508	31	4.596	0.000

plot(t1.rt.m1, horizontal = FALSE, comparisons = T)

t1.rt.m2 <- emmeans(t1.rt.aov1, pairwise ~ cCong | cnRep) # adjust="bon" , "tukey"
t1.rt.m2$contrasts %>% kable(digits=3)

contrast	cnRep	estimate	SE	df	t.ratio	p.value
v - nv	X1	-1.966	6.290	31	-0.313	0.757
v - nv	X0	-1.165	7.615	31	-0.153	0.879

plot(t1.rt.m2, horizontal = FALSE, comparisons = T)

2.1.2 RT analysis: Priming

t3.0 <- t1 %>% filter(Corr == 1) 

t3.1 <- t3.0 %>% filter(RT > 200 & RT < 10000) %>%
  group_by(SN) %>% # grouping by participants
  nest() %>%
  mutate(lbound = map(data, ~mean(.$RT)-3*sd(.$RT)),
         ubound = map(data, ~mean(.$RT)+3*sd(.$RT))) %>% # make new data (3sd cut)
  unnest(c(lbound, ubound))%>%
  unnest(data) %>%
  mutate(Outlier = (RT < lbound)|(RT > ubound)) %>% # set outlier
  filter(Outlier == FALSE) %>% # filtering outlier
  ungroup()

# 100-100*(nrow(t3.1)/nrow(t3.0))

t3.1 <- t3.1 %>% filter(cnRep!=0, cnRep!=1, cnRep!=2)
glimpse(t3.1, width = 50)
## Rows: 10,873
## Columns: 18
## $ SN      <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ Trial   <int> 25, 26, 27, 29, 30, 31, 32, 34, …
## $ Block   <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ bTrial  <int> 25, 26, 27, 29, 30, 31, 32, 34, …
## $ cCong   <fct> nv, nv, v, v, v, nv, nv, nv, nv,…
## $ cnRep   <fct> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4,…
## $ cnPair  <fct> 1, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2,…
## $ cIdt    <fct> 4, 4, 2, 5, 5, 3, 3, 6, 6, 1, 4,…
## $ cCat    <fct> 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 1,…
## $ Im_nExm <int> 3, 1, 3, 3, 1, 3, 1, 3, 1, 1, 1,…
## $ Im_Idx  <fct> 209, 111, 47, 181, 72, 157, 127,…
## $ Im_Name <fct> 209_3_out.jpg, 111_1_in.jpg, 047…
## $ Resp    <int> 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 1,…
## $ RT      <dbl> 589.8, 546.3, 935.2, 591.1, 525.…
## $ Corr    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ lbound  <dbl> 199.812, 199.812, 199.812, 199.8…
## $ ubound  <dbl> 996.88, 996.88, 996.88, 996.88, …
## $ Outlier <lgl> FALSE, FALSE, FALSE, FALSE, FALS…

t3.2 <- t3.1 %>% filter(cnPair == 2) %>% 
  dplyr::select(SN, cCong, cnRep, cIdt, Im_Name, RT)

t3.3 <- t3.2 %>% group_by(SN) %>% 
  spread(key = cnRep, value = RT) %>% ungroup()

t3.4 <- t3.3 %>% group_by(SN, cCong, cIdt) %>% mutate(prime = `4` - `3`) %>% ungroup()
t3.4 <- na.omit(t3.4)
glimpse(t3.4)
## Rows: 3,505
## Columns: 7
## $ SN      <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ cCong   <fct> v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v,…
## $ cIdt    <fct> 1, 5, 8, 11, 12, 13, 14, 15, 16, 21, 22, 26, 27, 30, 31, 33, 3…
## $ Im_Name <fct> 193_1_out.jpg, 072_1_in.jpg, 124_1_in.jpg, 150_1_out.jpg, 212_…
## $ `3`     <dbl> 520.2, 525.2, 755.2, 722.3, 487.2, 515.7, 552.4, 518.7, 535.0,…
## $ `4`     <dbl> 471.9, 560.4, 588.6, 549.8, 568.8, 513.3, 536.6, 551.9, 522.2,…
## $ prime   <dbl> -48.3, 35.2, -166.6, -172.5, 81.6, -2.4, -15.8, 33.2, -12.8, -…

# subject-level, long format
t3rtL <- t3.4 %>% group_by(SN, cCong) %>%
  summarise(prime = mean(prime)) %>%
  ungroup()
## `summarise()` has grouped output by 'SN'. You can override using the `.groups` argument.
#t3rtL %>% kable(digits=2)

# subject-level, wide format (SN/Btw/Block)
t3rtW <- t3rtL %>% spread(key=cCong, value = prime)
#t3rtW %>% kable(digits=2)

# summary table: grand mean (eyerep/locrep)
t3rtG <- t3rtL %>% group_by(cCong) %>%
  summarise(prime.m = mean(prime), prime.sd = sd(prime)) %>%
  ungroup()
t3rtG$prime.se <- Rmisc::summarySEwithin(data = t3rtL, measurevar = "prime", 
                                         idvar = "SN", withinvars = c("cCong"))$se
t3rtG$prime.ci <- Rmisc::summarySEwithin(data = t3rtL, measurevar = "prime", 
                                         idvar = "SN", withinvars = c("cCong"))$ci
t3rtG <- t3rtG %>% 
  mutate(lower.ci = prime.m-prime.ci,
         upper.ci = prime.m+prime.ci)
t3rtG %>% kable(digits=2)

cCong	prime.m	prime.sd	prime.se	prime.ci	lower.ci	upper.ci
v	10.68	42.24	8.41	17.16	-6.48	27.84
nv	-16.03	41.28	8.41	17.16	-33.19	1.13

# plot
p1.L.rt3.g <- t3rtG
p1.L.rt3.g$Cong <- factor(p1.L.rt3.g$cCong, labels=c("Violation", "Non-violation"))
p1.L.rt3.l <- t3rtL
p1.L.rt3.l$Cong <- factor(p1.L.rt3.l$cCong, labels=c("Violation", "Non-violation"))
p1.L.rt3.w <- p1.L.rt3.l %>% dplyr::select(SN, Cong, prime) %>% spread(key =Cong, value=prime)

p1.L.rt3.g.plot3 <- ggplot() +
  geom_hline(yintercept=0, linetype='dashed', color='gray60', alpha =1, size=1) +
  geom_point(data=p1.L.rt3.l, aes(x=Cong, y=prime), size = 1, position = position_dodge(.8),  color="gray80") +
  geom_segment(data=p1.L.rt3.w, inherit.aes = FALSE,
               aes(x=1, y=`Violation`,
                   xend=2, yend=`Non-violation`),
               color="gray80") +
  geom_pointrange(data=p1.L.rt3.g, aes(x = Cong, y=prime.m, ymin = prime.m-prime.ci, ymax = prime.m+prime.ci, color = Cong),
                  position = position_dodge(0.80), size = 1.5, show.legend = FALSE) +
  scale_color_manual(values = c("#ED7D31", "#70AD47")) + 
  # scale_shape_manual(values = c("red", "black"), labels = c("Experimental", "Control")) +
  coord_cartesian(ylim = c(-150,+150), clip = "on") +
  # scale_y_continuous(breaks = seq(0,80)) +
  scale_x_discrete(labels = c("Violation", "Non-violation")) +
  labs(x = "Type", y = "4th - 3rd (ms)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        legend.position = c(0.65, 0.85),
        legend.title = element_blank(),
        # aspect.ratio = 0.6,
        legend.background = element_blank(),
        plot.margin = margin(1, 0.3, 0.3, 0.5, "cm"),
        legend.key = element_blank())

p1.L.rt3.g.plot3

# p1.rt.plot2 <- ggarrange(p1.L.rt.g.plot1, p1.L.rt.g.plot2,p1.L.rt3.g.plot3, ncol = 3,
#                          labels=c("A) Categorization of B", "B) Categorization fo B/B`", "C) Differentiation Effect"))
# p1.rt.plot2

p1.L.rt3.g.plot3.sn <- ggplot() +
  geom_hline(yintercept=0,  color='gray60', alpha =1, size=1) +
  stat_summary(data=p1.L.rt3.l, aes(x=SN, y=prime, fill=Cong), 
               fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, 
               width = 0.8,size = 1.02) + # , show.legend = FALSE  colour="black", 
  
  scale_fill_manual(values = c("#ED7D31", "#70AD47"),
                    labels = c("Violation", "Non-violation")) + 
  # scale_shape_manual(values = c("red", "black"), labels = c("Experimental", "Control")) +
  coord_cartesian(ylim = c(-120,+120), clip = "on") +
  # scale_y_continuous(breaks = seq(0,80)) +
  # scale_x_discrete(labels = c("Violation", "Non-violation")) +
  labs(x = "Type", y = "4th - 3rd (ms)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        legend.position = 'top', # c(0.65, 0.85),
        legend.title = element_blank(),
        # aspect.ratio = 0.6,
        legend.background = element_blank(),
        plot.margin = margin(1, 0.3, 0.3, 0.5, "cm"),
        legend.key = element_blank())

p1.L.rt3.g.plot3.sn

Paired T-Test 1 - violation 여부에 따른 priming 약화 효과 (differentiation effect)

t3rtL.1 <- t3rtL
t3.rt.ttest <- t.test(prime ~ cCong, data=t3rtL.1,
                      alternative = c("two.sided"),
                      paired=T, var.equal=F,
                      conf.level=.95)
t3.rt.ttest

## 
##  Paired t-test
## 
## data:  prime by cCong
## t = 2.2441, df = 31, p-value = 0.03211
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.43423 50.97299
## sample estimates:
## mean of the differences 
##                26.70361

t3.rt.rmanova <- aov_ez(id="SN", dv = "prime", data = t3rtL.1, within = c("cCong"))
nice(t3.rt.rmanova, es ="pes")

## Anova Table (Type 3 tests)
## 
## Response: prime
##   Effect    df     MSE      F  pes p.value
## 1  cCong 1, 31 2265.60 5.04 * .140    .032
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

t3.rt.m1 <- emmeans(t3.rt.rmanova, pairwise ~ cCong, adjust="tukey") # adjust="bon" , "tukey"
t3.rt.m1$contrasts %>% kable(digits=3)

contrast	estimate	SE	df	t.ratio	p.value
v - nv	26.704	11.9	31	2.244	0.032

plot(t3.rt.m1, horizontal = FALSE, comparisons = T)

2.1.3 RT analysis: Priming with B’ RT

t4.0 <- t1 %>% filter(Corr == 1)

t4.1 <- t4.0 %>% filter(RT > 200 & RT < 10000) %>%
  group_by(SN) %>% # grouping by participants
  nest() %>%
  mutate(lbound = map(data, ~mean(.$RT)-3*sd(.$RT)),
         ubound = map(data, ~mean(.$RT)+3*sd(.$RT))) %>% # make new data (3sd cut)
  unnest(c(lbound, ubound))%>%
  unnest(data) %>%
  mutate(Outlier = (RT < lbound)|(RT > ubound)) %>% # set outlier
  filter(Outlier == FALSE) %>% # filtering outlier
  ungroup()

100-100*(nrow(t4.1)/nrow(t4.0))
## [1] 1.812514

t4.1 <- t4.1 %>% filter(cnRep!=1, cnRep!=2)
# t4.1$cnRep <- as.double(t4.1$cnRep)
t4.2 <- t4.1 %>% filter(cnPair == 2) %>% 
  dplyr::select(SN, cCong, cnRep, cIdt, Im_Idx, RT)

t4.3 <- t4.2 %>% group_by(SN) %>% 
  spread(key = cnRep, value = RT) %>% ungroup()


t4.4 <- t4.3 %>% group_by(SN, cIdt) %>% mutate(prime = `4` - `3`) %>% ungroup()
t4.4 <- na.omit(t4.4)
t4.4 <- t4.4 %>% dplyr::select (SN, cCong, cIdt, Im_Idx, `0`, prime)
glimpse(t4.4)
## Rows: 3,345
## Columns: 6
## $ SN     <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ cCong  <fct> v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, …
## $ cIdt   <fct> 1, 5, 11, 12, 13, 14, 16, 21, 22, 26, 27, 30, 31, 33, 34, 35, 3…
## $ Im_Idx <fct> 193, 72, 150, 212, 234, 220, 122, 149, 109, 180, 207, 95, 80, 7…
## $ `0`    <dbl> 663.1, 631.0, 596.8, 548.1, 534.3, 752.2, 516.5, 685.6, 572.1, …
## $ prime  <dbl> -48.3, 35.2, -172.5, 81.6, -2.4, -15.8, -12.8, -92.9, -21.5, -6…

# subject-level, long format
t4rtL <- t4.4 %>% group_by(SN, cCong) %>%
  summarise(prime = mean(prime), bpr = mean(`0`)) %>%
  ungroup()
## `summarise()` has grouped output by 'SN'. You can override using the `.groups` argument.
#t4rtL %>% kable(digits=2)

# summary table: grand mean (eyerep/locrep)
t4rtG <- t4rtL %>% group_by(cCong) %>%
  summarise(prime.m = mean(prime), prime.sd = sd(prime), bpr.m = mean(bpr), brp.sd =sd(bpr)) %>%
  ungroup()
t3rtG %>% kable(digits=2)

cCong	prime.m	prime.sd	prime.se	prime.ci	lower.ci	upper.ci
v	10.68	42.24	8.41	17.16	-6.48	27.84
nv	-16.03	41.28	8.41	17.16	-33.19	1.13

# correlation by Subject
t4rtL.1 <- t4rtL %>% filter(cCong == 'v')
cor.test(formula = ~ prime + bpr, data = t4rtL.1,
         method = "pearson", alternative = "two.sided")

## 
##  Pearson's product-moment correlation
## 
## data:  prime and bpr
## t = 1.0021, df = 30, p-value = 0.3243
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1800267  0.4974426
## sample estimates:
##       cor 
## 0.1799644

t4rtL.2 <- t4rtL %>% filter(cCong == 'nv')
cor.test(formula = ~ prime + bpr, data = t4rtL.2,
         method = "pearson", alternative = "two.sided")

## 
##  Pearson's product-moment correlation
## 
## data:  prime and bpr
## t = 0.7508, df = 30, p-value = 0.4586
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2234696  0.4625946
## sample estimates:
##       cor 
## 0.1358068

p1.rt.corr.plot1 <- ggarrange(p1.L.rt.corr.plot1, p1.L.rt.corr.plot2, ncol = 2,
                         labels=c("A) Violation Condition", "B) Non-violation Condition"))

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

p1.rt.corr.plot1

2.1.4 RT analysis: correlation, aggregation

t4.4.1 <- t4.4 %>% filter(cCong =="v")
cor.test(formula = ~ prime + `0`, data = t4.4.1,
         method = "pearson", alternative = "two.sided")

## 
##  Pearson's product-moment correlation
## 
## data:  prime and 0
## t = 0.63219, df = 1646, p-value = 0.5273
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.03273072  0.06381923
## sample estimates:
##        cor 
## 0.01558057

ggplot(t4.4.1, aes(x=`0`, y=prime)) +
  geom_point(size = 4) +
  geom_smooth(method=lm) +
  labs(x = "b` RT", 
       y = "4th - 3rd B RT") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1)

## `geom_smooth()` using formula 'y ~ x'

t4.4.2 <- t4.4 %>% filter(cCong =="nv")
cor.test(formula = ~ prime + `0`, data = t4.4.2,
         method = "pearson", alternative = "two.sided")

## 
##  Pearson's product-moment correlation
## 
## data:  prime and 0
## t = 2.0711, df = 1695, p-value = 0.0385
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.002663411 0.097592241
## sample estimates:
##       cor 
## 0.0502413

ggplot(t4.4.2, aes(x=`0`, y=prime)) +
  geom_point(size = 4) +
  geom_smooth(method=lm) +
  labs(x = "b` RT", 
       y = "4th - 3rd B RT") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1)

## `geom_smooth()` using formula 'y ~ x'

2.1.5 Accuracy analysis

t1acc.sn <- t1 %>% filter(cnPair == 2) %>% group_by(SN, cCong) %>%
  dplyr::summarise(Acc = mean(Corr)*100) %>%
  ungroup() # %>% spread(key=cCong,value=Acc)
## `summarise()` has grouped output by 'SN'. You can override using the `.groups` argument.
#

ggplot(data=t1acc.sn, aes(x = SN, y = Acc)) + 
  ggtitle("mean Accuracy") +
  xlab("Participant") + ylab("Acc") + theme_bw() +
  theme(text=element_text(size=14)) +
  theme(legend.title = element_text(face = 1,size = 15)) +
  geom_boxplot()

t1accL <- t1 %>% filter(cnPair == 2) %>% 
  group_by(SN, cCong, cnRep) %>%
  dplyr::summarise(Acc = mean(Corr)*100) %>%
  ungroup()
## `summarise()` has grouped output by 'SN', 'cCong'. You can override using the `.groups` argument.
# t1accL %>% kable(digits=2)

# subject-level, wide format (SN/Btw/Block)
t1accW1 <- t1accL %>% spread(key=cCong, value = Acc)
# t1accW1 %>% kable(digits=2)

t1accW2 <- t1accL %>% spread(key=cnRep, value = Acc)
# t1accW2 %>% kable(digits=2)

# summary table: grand mean (eyerep/locrep)
t1accG <- t1accL %>% group_by(cCong, cnRep) %>%
  summarise(Acc.m = mean(Acc), Acc.sd = sd(Acc)) %>%
  ungroup()
## `summarise()` has grouped output by 'cCong'. You can override using the `.groups` argument.
t1accG$Acc.se <- Rmisc::summarySEwithin(data = t1accL, measurevar = "Acc", 
                                        idvar = "SN", withinvars = c("cCong", "cnRep"))$se
t1accG$Acc.ci <- Rmisc::summarySEwithin(data = t1accL, measurevar = "Acc", 
                                        idvar = "SN", withinvars = c("cCong", "cnRep"))$ci
t1accG <- t1accG %>% 
  mutate(lower.ci = Acc.m-Acc.ci,
         upper.ci = Acc.m+Acc.ci)
t1accG %>% 
  dplyr::select(cCong, cnRep, Acc.m) %>% 
  spread(key=cnRep, value= Acc.m) %>% kable(digits=2)

cCong	1	2	3	4	0
v	93.85	95.16	95.42	95.21	95.57
nv	94.84	95.42	96.35	96.51	96.30

<br

p1.L.acc.g <- t1accG %>% filter(cnRep != 0)
p1.L.acc.g$Cong <- factor(p1.L.acc.g$cCong, labels=c("Violation", "Non-violation"))
p1.L.acc.g.plot1 <- ggplot(data=p1.L.acc.g , 
                           aes(x=cnRep, y=Acc.m, ymin=Acc.m-Acc.ci, ymax=Acc.m+Acc.ci, color=Cong, shape=Cong)) +
  geom_point(size = 5, position = position_dodge(.3)) +
  geom_errorbar(width = .2, position = position_dodge(.3)) +
  geom_line(aes(group = Cong), position = position_dodge(.3), size=1) + 
  scale_color_manual(values = c("#ED7D31", "#70AD47")) + # c("red", "black")) +
  # scale_shape_manual(values = c("red", "black"), labels = c("Experimental", "Control")) +
  coord_cartesian(ylim = c(80,110), clip = "on") +
  # scale_y_continuous(breaks = seq(0,80)) +
  scale_x_discrete(labels = c("1", "2", "3", "4")) +
  labs(x = "Repetition", y = "Accuracy (%)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        legend.position = c(0.65, 0.85),
        legend.title = element_blank(),
        # aspect.ratio = 0.6,
        legend.background = element_blank(),
        plot.margin = margin(1, 0.3, 0.3, 0.5, "cm"),
        legend.key = element_blank())

p1.L.acc.g <- t1accG %>% filter(cnRep != 2, cnRep != 3, cnRep != 4)
p1.L.acc.g$Cong <- factor(p1.L.acc.g$cCong, labels=c("Violation", "Non-violation"))
p1.L.acc.g.plot2 <- ggplot(data=p1.L.acc.g , 
                           aes(x=cnRep, y=Acc.m, ymin=Acc.m-Acc.ci, ymax=Acc.m+Acc.ci, color=Cong, shape=Cong)) +
  geom_point(size = 5, position = position_dodge(.3)) +
  geom_errorbar(width = .2, position = position_dodge(.3)) +
  geom_line(aes(group = Cong), position = position_dodge(.3), size=1) + 
  scale_color_manual(values = c("#ED7D31", "#70AD47")) + # c("red", "black")) +
  # scale_shape_manual(values = c("red", "black"), labels = c("Experimental", "Control")) +
  coord_cartesian(ylim = c(80,110), clip = "on") +
  # scale_y_continuous(breaks = seq(0,80)) +
  scale_x_discrete(labels = c("1st B", "B-Prime")) +
  labs(x = "Type", y = "Accuracy (%)") +
  theme_bw(base_size = 18) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        legend.position = c(0.65, 0.85),
        legend.title = element_blank(),
        # aspect.ratio = 0.6,
        legend.background = element_blank(),
        plot.margin = margin(1, 0.3, 0.3, 0.5, "cm"),
        legend.key = element_blank())

p1.acc.plot <- ggarrange(p1.L.acc.g.plot1, p1.L.acc.g.plot2, ncol = 2,
                         labels=c("A) Categorization of B", "B) Categorization fo B/B`"))
p1.acc.plot

RM ANOVA

t1accL.1 <- t1accL # %>% filter(SN!=4) 
t1.acc.aov1 <- aov_ez(id="SN", dv = "Acc", data = t1accL.1, within = c("cCong", "cnRep"))
# test_sphericity(t2.acc.aov1)
#summary(t1.acc.aov1)
#anova(t1.acc.aov1, es = "pes")
nice(t1.acc.aov1, es="pes") %>% kable(digits=3)

Effect	df	MSE	F	pes	p.value
cCong	1, 31	14.98	3.80 +	.109	.060
cnRep	3.09, 95.74	7.73	4.91 **	.137	.003
cCong:cnRep	3.16, 97.87	6.57	0.46	.015	.722

RM ANOVA 1 - post-hoc

t1.acc.m1 <- emmeans(t1.acc.aov1, pairwise ~ cnRep | cCong) # adjust="bon" , "tukey"
t1.acc.m1$contrasts %>% kable(digits=3)

contrast	cCong	estimate	SE	df	t.ratio	p.value
X1 - X2	v	-1.302	0.634	31	-2.054	0.265
X1 - X3	v	-1.562	0.617	31	-2.533	0.109
X1 - X4	v	-1.354	0.582	31	-2.328	0.163
X1 - X0	v	-1.719	0.615	31	-2.796	0.062
X2 - X3	v	-0.260	0.627	31	-0.416	0.993
X2 - X4	v	-0.052	0.663	31	-0.079	1.000
X2 - X0	v	-0.417	0.594	31	-0.701	0.955
X3 - X4	v	0.208	0.593	31	0.351	0.997
X3 - X0	v	-0.156	0.542	31	-0.289	0.998
X4 - X0	v	-0.365	0.528	31	-0.691	0.957
X1 - X2	nv	-0.573	0.659	31	-0.869	0.906
X1 - X3	nv	-1.510	0.542	31	-2.789	0.063
X1 - X4	nv	-1.667	0.698	31	-2.388	0.146
X1 - X0	nv	-1.458	0.701	31	-2.080	0.254
X2 - X3	nv	-0.937	0.449	31	-2.090	0.250
X2 - X4	nv	-1.094	0.583	31	-1.877	0.350
X2 - X0	nv	-0.885	0.467	31	-1.895	0.341
X3 - X4	nv	-0.156	0.493	31	-0.317	0.998
X3 - X0	nv	0.052	0.505	31	0.103	1.000
X4 - X0	nv	0.208	0.634	31	0.329	0.997

plot(t1.acc.m1, horizontal = FALSE, comparisons = T)

t1.acc.m2 <- emmeans(t1.acc.aov1, pairwise ~ cCong | cnRep) # adjust="bon" , "tukey"
t1.acc.m2$contrasts %>% kable(digits=3)

contrast	cnRep	estimate	SE	df	t.ratio	p.value
v - nv	X1	-0.990	0.681	31	-1.453	0.156
v - nv	X2	-0.260	0.508	31	-0.512	0.612
v - nv	X3	-0.937	0.639	31	-1.467	0.152
v - nv	X4	-1.302	0.829	31	-1.571	0.126
v - nv	X0	-0.729	0.643	31	-1.133	0.266

plot(t1.acc.m2, horizontal = FALSE, comparisons = T)

RM-ANOVA 2

t1accL.2 <- t1accL %>% filter(cnRep !=1, cnRep != 2, cnRep !=0)
t1.acc.aov1 <- aov_ez(id="SN", dv = "Acc", data = t1accL.2, within = c("cCong", "cnRep"))
# test_sphericity(t2.acc.aov1)
# summary(t1.acc.aov1)
# anova(t1.acc.aov1, es = "pes")
nice(t1.acc.aov1, es="pes") %>% kable(digits=3)

Effect	df	MSE	F	pes	p.value
cCong	1, 31	11.99	3.35 +	.097	.077
cnRep	1, 31	3.96	0.01	<.001	.941
cCong:cnRep	1, 31	5.54	0.19	.006	.664

RM ANOVA 2 - post-hoc

t1.acc.m1 <- emmeans(t1.acc.aov1, pairwise ~ cCong | cnRep) # adjust="bon" , "tukey"
t1.acc.m1$contrasts %>% kable(digits=3)

contrast	cnRep	estimate	SE	df	t.ratio	p.value
v - nv	X3	-0.937	0.639	31	-1.467	0.152
v - nv	X4	-1.302	0.829	31	-1.571	0.126

plot(t1.acc.m1, horizontal = FALSE, comparisons = T)

t1.acc.m1 <- emmeans(t1.acc.aov1, pairwise ~ cnRep) # adjust="bon" , "tukey"
t1.acc.m1$contrasts %>% kable(digits=3)

contrast	estimate	SE	df	t.ratio	p.value
X3 - X4	0.026	0.352	31	0.074	0.941

plot(t1.acc.m1, horizontal = FALSE, comparisons = T)

Paired T-Test

t1accL.3 <- t1accL %>% filter(cnRep == 0)
t1.acc.ttest <- t.test(Acc ~ cCong, data=t1accL.3,
                       alternative = c("two.sided"),
                       paired=T, var.equal=F,
                       conf.level=.95)
t1.acc.ttest 

## 
##  Paired t-test
## 
## data:  Acc by cCong
## t = -1.1331, df = 31, p-value = 0.2658
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.0415686  0.5832353
## sample estimates:
## mean of the differences 
##              -0.7291667

RM-ANOVA 3

t1accL.3 <- t1accL %>% filter(cnRep !=2, cnRep !=3, cnRep !=4)
t1.acc.aov1 <- aov_ez(id="SN", dv = "Acc", data = t1accL.3, within = c("cCong", "cnRep"))
# test_sphericity(t2.acc.aov1)
# summary(t1.acc.aov1)
# anova(t1.acc.aov1, es = "pes")
nice(t1.acc.aov1, es="pes") %>% kable(digits=3)

Effect	df	MSE	F	pes	p.value
cCong	1, 31	9.79	2.41	.072	.130
cnRep	1, 31	9.65	8.37 **	.213	.007
cCong:cnRep	1, 31	4.26	0.13	.004	.724

RM-ANOVA 3 - post hoc

### RM ANOVA 3 - post-hoc
t1.acc.m1 <- emmeans(t1.acc.aov1, pairwise ~ cnRep | cCong) # adjust="bon" , "tukey"
t1.acc.m1$contrasts %>% kable(digits=3)

contrast	cCong	estimate	SE	df	t.ratio	p.value
X1 - X0	v	-1.719	0.615	31	-2.796	0.009
X1 - X0	nv	-1.458	0.701	31	-2.080	0.046

plot(t1.acc.m1, horizontal = FALSE, comparisons = T)

t1.rt.m2 <- emmeans(t1.rt.aov1, pairwise ~ cCong | cnRep) # adjust="bon" , "tukey"
t1.rt.m2$contrasts %>% kable(digits=3)

contrast	cnRep	estimate	SE	df	t.ratio	p.value
v - nv	X1	-1.966	6.290	31	-0.313	0.757
v - nv	X0	-1.165	7.615	31	-0.153	0.879

plot(t1.rt.m2, horizontal = FALSE, comparisons = T)

---

2.2 Phase 2 - Recognition Test

t2 <- read.csv("data/reactMEM_p2_all.csv", header = T)

glimpse(t2, width = 70)
## Rows: 7,680
## Columns: 13
## $ SN       <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ Trial    <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, …
## $ cType    <int> 2, 1, 1, 2, 1, 2, 2, 1, 2, 2, 1, 2, 2, 1, 2, 1, 2, …
## $ cCong    <int> 999, 1, 2, 999, 1, 999, 999, 2, 999, 999, 2, 999, 9…
## $ cIdt     <int> 999, 117, 66, 999, 60, 999, 999, 95, 999, 999, 113,…
## $ cCat     <int> 2, 2, 1, 1, 2, 2, 2, 1, 2, 1, 2, 2, 1, 1, 2, 1, 2, …
## $ xIm_nExm <int> 999, 2, 2, 999, 2, 999, 999, 2, 999, 999, 2, 999, 9…
## $ xIm_Idx  <int> 241, 227, 71, 301, 191, 341, 303, 27, 252, 282, 143…
## $ dIm_Name <chr> "241_1_out.jpg", "227_2_out.jpg", "071_2_in.jpg", "…
## $ Resp     <int> 1, 4, 2, 2, 2, 1, 1, 4, 1, 2, 4, 2, 1, 3, 2, 2, 1, …
## $ RT       <dbl> 2.8079, 2.2212, 2.0978, 1.2932, 1.6856, 1.5962, 1.3…
## $ Corr     <int> 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, …
## $ Conf     <int> 2, 2, 1, 1, 1, 2, 2, 2, 2, 1, 2, 1, 2, 1, 1, 1, 2, …

# ------- Variables ------- #
# $ SN       <int> 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99…
# $ Trial    <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,…
# $ cType    <int> 1, 1, 1, 2, 2, 1, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1,…
# $ cCong    <int> 1, 2, 1, 999, 999, 1, 999, 999, 2, 1, 1, 1, 999, 9…
# $ cIdt     <int> 3, 105, 43, 999, 999, 89, 999, 999, 103, 81, 23, 2…
# $ cCat     <int> 2, 1, 2, 2, 2, 1, 2, 1, 2, 1, 1, 1, 1, 1, 2, 2, 1,…
# $ xIm_nExm <int> 2, 2, 2, 999, 999, 2, 999, 999, 2, 2, 2, 2, 999, 9…
# $ xIm_Idx  <int> 199, 130, 197, 271, 259, 4, 278, 274, 192, 73, 62,…
# $ dIm_Name <chr> " 199_2_out.jpg", " 130_2_in.jpg", " 197_2_out.jpg…
# $ Resp     <int> 2, 4, 3, 1, 1, 4, 4, 1, 2, 4, 2, 1, 1, 1, 1, 3, 4,…
# $ RT       <dbl> 62.2271, 2.4708, 2.2955, 5.4393, 1.4200, 1.1727, 1…
# $ Corr     <int> 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1,…
# $ Conf     <int> 1, 2, 1, 2, 2, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 1, 2,…

length(unique(t2$SN))
## [1] 32

# check number of trials for each condition/SN
table(t2$SN)
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
## 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 
##  21  22  23  24  25  26  27  28  29  30  31  32 
## 240 240 240 240 240 240 240 240 240 240 240 240
table(t2$cType, t2$SN) 
##    
##       1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19
##   1 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120
##   2 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120
##    
##      20  21  22  23  24  25  26  27  28  29  30  31  32
##   1 120 120 120 120 120 120 120 120 120 120 120 120 120
##   2 120 120 120 120 120 120 120 120 120 120 120 120 120
table(t2$cCong, t2$SN) 
##      
##         1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
##   1    60  60  60  60  60  60  60  60  60  60  60  60  60  60  60  60  60  60
##   2    60  60  60  60  60  60  60  60  60  60  60  60  60  60  60  60  60  60
##   999 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120
##      
##        19  20  21  22  23  24  25  26  27  28  29  30  31  32
##   1    60  60  60  60  60  60  60  60  60  60  60  60  60  60
##   2    60  60  60  60  60  60  60  60  60  60  60  60  60  60
##   999 120 120 120 120 120 120 120 120 120 120 120 120 120 120
# table(t2$cCat, t2$SN) 

# change class of main factors: double to factor
t2$SN = factor(t2$SN)
t2$cType = factor(t2$cType, levels=c(1,2), labels=c("old","new"))
t2$cCong = factor(t2$cCong, levels=c(1,2,999), labels=c("v","nv","new")) 
t2$cIdt = factor(t2$cIdt)
# t2$cCat = factor(t2$cCat)
t2$dIm_Name = factor(t2$dIm_Name)
t2$Resp <- factor(t2$Resp);
t2$mResp <- factor(t2$Resp, levels=c(1,2,3,4), labels=c("new", "new", "old", "old"))
t2$RT <- t2$RT*1000;
t2$Corr <- as.numeric(t2$Corr==1)
t2$Conf <- as.double(t2$Conf)
# t2$Corr[t2$Conf==2] <- 1
# t2$Corr[t2$cType=="old" & t2$mResp=="old" &  t2$Conf==2] <- 1
# t2$Corr[t2$cType=="old" & t2$mResp=="old" & t2$Conf==1] <- 0

t2$p_resp = factor(t2$cType, levels=c('old','new'), labels=c("p.old","p.new"))
t2$r_resp = factor(t2$mResp, levels=c('old','new'), labels=c("r.old", "r.new"))

glimpse(t2)
## Rows: 7,680
## Columns: 16
## $ SN       <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ Trial    <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
## $ cType    <fct> new, old, old, new, old, new, new, old, new, new, old, new, n…
## $ cCong    <fct> new, v, nv, new, v, new, new, nv, new, new, nv, new, new, nv,…
## $ cIdt     <fct> 999, 117, 66, 999, 60, 999, 999, 95, 999, 999, 113, 999, 999,…
## $ cCat     <int> 2, 2, 1, 1, 2, 2, 2, 1, 2, 1, 2, 2, 1, 1, 2, 1, 2, 2, 2, 2, 2…
## $ xIm_nExm <int> 999, 2, 2, 999, 2, 999, 999, 2, 999, 999, 2, 999, 999, 2, 999…
## $ xIm_Idx  <int> 241, 227, 71, 301, 191, 341, 303, 27, 252, 282, 143, 280, 245…
## $ dIm_Name <fct> 241_1_out.jpg, 227_2_out.jpg, 071_2_in.jpg, 301_1_in.jpg, 191…
## $ Resp     <fct> 1, 4, 2, 2, 2, 1, 1, 4, 1, 2, 4, 2, 1, 3, 2, 2, 1, 2, 4, 2, 2…
## $ RT       <dbl> 2807.9, 2221.2, 2097.8, 1293.2, 1685.6, 1596.2, 1394.6, 1831.…
## $ Corr     <dbl> 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1…
## $ Conf     <dbl> 2, 2, 1, 1, 1, 2, 2, 2, 2, 1, 2, 1, 2, 1, 1, 1, 2, 1, 2, 1, 1…
## $ mResp    <fct> new, old, new, new, new, new, new, old, new, new, old, new, n…
## $ p_resp   <fct> p.new, p.old, p.old, p.new, p.old, p.new, p.new, p.old, p.new…
## $ r_resp   <fct> r.new, r.old, r.new, r.new, r.new, r.new, r.new, r.old, r.new…

2.2.1 Accuracy Analysis: Raw Data

t2acc.sn <- t2 %>% group_by(SN, cType, cCong) %>%
  dplyr::summarise(Acc = mean(Corr)*100) %>%
  ungroup() %>% spread(key=cType,value=Acc)

## `summarise()` has grouped output by 'SN', 'cType'. You can override using the `.groups` argument.

# t2acc.sn %>% kable(digits=2)

t2acc.snL <- t2 %>% 
  filter(cType=='old') %>%
  group_by(SN, cCong) %>%
  dplyr::summarise(Acc = mean(Corr)*100) %>%
  ungroup()

## `summarise()` has grouped output by 'SN'. You can override using the `.groups` argument.

#t2acc.snL %>% kable(digits=2)

ggplot(data=t2acc.snL, aes(x = SN, y = Acc)) + 
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, 
               width = 0.8,size = 1.02) + # , show.legend = FALSE  colour="black",
  geom_hline(yintercept=50, linetype='dashed', color='darkred', alpha =1, size=1) +
  coord_cartesian(ylim = c(0, 100), clip = "on") +
  ggtitle("mean Accuracy (studied)") +
  xlab("Participants") + ylab("Acc") + theme_bw() +
  theme(text=element_text(size=14)) +
  theme(legend.title = element_text(face = 1,size = 15))

# subject-level, long format (SN/Btw/Block)
t2accL <- t2 %>% group_by(SN, cType, cCong) %>%
  dplyr::summarise(Acc = mean(Corr)*100) %>%
  ungroup()
## `summarise()` has grouped output by 'SN', 'cType'. You can override using the `.groups` argument.
# t2accL %>% kable(digits=2)

# subject-level, wide format (SN/Btw/Block)
t2accW1 <- t2accL %>% spread(key=cCong, value = Acc)
# t2accW1 %>% kable(digits=2)

t2accW2 <- t2accL %>% spread(key=cType, value = Acc)
# t2accW2 %>% kable(digits=2)

# summary table: grand mean (eyerep/locrep)
t2accG <- t2accL %>% group_by(cType, cCong) %>%
  summarise(Acc.m = mean(Acc), Acc.sd = sd(Acc)) %>%
  ungroup()
## `summarise()` has grouped output by 'cType'. You can override using the `.groups` argument.
t2accG$Acc.se <- Rmisc::summarySEwithin(data = t2accL, measurevar = "Acc", 
                                        idvar = "SN", withinvars = c("cType", "cCong"))$se
t2accG$Acc.ci <- Rmisc::summarySEwithin(data = t2accL, measurevar = "Acc", 
                                        idvar = "SN", withinvars = c("cType", "cCong"))$ci
t2accG <- t2accG %>% 
  mutate(lower.ci = Acc.m-Acc.ci,
         upper.ci = Acc.m+Acc.ci)
# for between-subject design. check help(summarySE)
t2accG %>% kable(digits=2)

cType	cCong	Acc.m	Acc.sd	Acc.se	Acc.ci	lower.ci	upper.ci
old	v	56.15	15.71	1.41	2.88	53.26	59.03
old	nv	52.40	16.34	1.51	3.09	49.31	55.48
new	new	94.45	6.25	2.51	5.12	89.34	99.57

t2.plot1 <- ggplot(data=t2accL, aes(x=cType, y=Acc, fill=cCong)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, 
               width = 0.8,size = 1.02) + # , show.legend = FALSE  colour="black", 
  geom_point(data=t2accL, aes(x=cType, y=Acc, fill=cCong), 
             position=position_dodge(width = 0.8), 
             size=2, show.legend = FALSE, color="gray60") + # , 
  geom_segment(data=filter(t2accW1, cType == "old"), inherit.aes = FALSE,
               aes(x=0.8, y=filter(t2accW1, cType == "old")$v,
                   xend=1.2, yend=filter(t2accW1, cType == "old")$nv),
               color="gray60") +
  geom_segment(data=filter(t2accW1, cType == "new"), inherit.aes = FALSE,
               aes(x=1.8, y=filter(t2accW1, cType == "new")$v,
                   xend=2.2, yend=filter(t2accW1, cType == "new")$nv),
               color="gray60") +
  geom_pointrange(data=t2accG, aes(x = cType, y=Acc.m, ymin = Acc.m-Acc.ci, ymax = Acc.m+Acc.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("Studied","Lure")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4", "#235796"),# c("#6798D4", "#D87B49", "#4FA77C"), 
                    labels = c("Violation", "Non-violation", "Lure")) +
  # scale_fill_manual(values = c("#feb24c", "#91bfdb"), #6798D4 #D87B49 #4FA77C
                    # labels = c("Violation", "Non-violation")) +
  coord_cartesian(ylim = c(0, 120), clip = "on") +
  scale_y_continuous(breaks = c(0,25,50,75,100)) +
  labs(x = "Stimuli Type", y = "Accuracy (%)", fill ="Violation Condition") +
  # ggtitle("Association Memory Accuracy") +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"), 
        # plot.title = element_text(hjust = 0.5),
        legend.position="top")
        # legend.position=c(0.8, 0.90))
t2.plot1

t2accL.sn <- t2accL %>% filter(cCong!="new")
t2.plot1.sn <- ggplot(data=t2accL.sn, aes(x=SN, y=Acc, fill=cCong)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, 
               width = 0.8,size = 1.02) + # , show.legend = FALSE  colour="black", 
  # scale_x_discrete(labels=c("Studied","Lure")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),# c("#6798D4", "#D87B49", "#4FA77C"), 
                    labels = c("Violation", "Non-violation")) +
  # scale_fill_manual(values = c("#feb24c", "#91bfdb"), #6798D4 #D87B49 #4FA77C
  # labels = c("Violation", "Non-violation")) +
  coord_cartesian(ylim = c(0, 120), clip = "on") +
  scale_y_continuous(breaks = c(0,25,50,75,100)) +
  labs(x = "Subject Number", y = "Accuracy (%)", fill ="Violation Condition") +
  # ggtitle("Association Memory Accuracy") +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"), 
        # plot.title = element_text(hjust = 0.5),
        legend.position="top")
# legend.position=c(0.8, 0.90))
t2.plot1.sn

2.2.2 Accuracy Analysis: Signal Detection Measure

t2sdtL.1.1 <- t2 %>% filter(cType=="old") %>% mutate( ) %>% group_by(SN, cCong, p_resp, r_resp) %>%
  summarise(n = sum(Corr==1 | Corr==0)/60) %>%
  unite(p_resp.r_resp, c(p_resp,r_resp), sep = "_", remove = TRUE) %>% ungroup() %>%
  spread(p_resp.r_resp, n)
## `summarise()` has grouped output by 'SN', 'cCong', 'p_resp'. You can override using the `.groups` argument.
t2sdtL.1.1 <- plyr::rename(t2sdtL.1.1, c("p.old_r.new" = "Miss", "p.old_r.old" = "Hit"))
t2sdtL.1.1 <- t2sdtL.1.1[,c(1, 2, 4, 3)]
# t2sdtL.1.1 %>% kable(digits=2)

t2sdtL.1.2 <- t2 %>% filter(cType=="new") %>% mutate( ) %>% group_by(SN, cCong, p_resp, r_resp) %>%
  # summarise(n = sum(Corr==1 | Corr==0)/60) %>%
  summarise(n = sum(Corr==1 | Corr==0)/120) %>%
  unite(p_resp.r_resp, c(p_resp,r_resp), sep = "_", remove = TRUE) %>% ungroup() %>%
  spread(p_resp.r_resp, n)
## `summarise()` has grouped output by 'SN', 'cCong', 'p_resp'. You can override using the `.groups` argument.
t2sdtL.1.2 <- plyr::rename(t2sdtL.1.2, c("p.new_r.new" = "CR", "p.new_r.old" = "FA"))
t2sdtL.1.2 <- t2sdtL.1.2[,c(1, 2, 4, 3)]
# t2sdtL.1.2 %>% kable(digits=2)


t2sdtL.3 <- merge(t2sdtL.1.1, t2sdtL.1.2, by="SN")
t2sdtL.3$cCong <- t2sdtL.3$cCong.x
t2sdtL.3[is.na(t2sdtL.3)] <-0
t2sdtL.3[,3] <- t2sdtL.3[,3] + 0.05
t2sdtL.3[,6] <- t2sdtL.3[,6] + 0.05
t2sdtL.3 <- t2sdtL.3 %>% dplyr::select(SN, cCong, Hit, Miss, FA, CR) %>%
  group_by(SN, cCong) %>% dplyr::mutate(mPerf = (Hit - FA)) %>% ungroup()
# t2sdtL.3 %>% kable(digits=2)

indices <- psycho::dprime(t2sdtL.3$Hit, t2sdtL.3$FA, n_targets=1, n_distractors = 1, adjusted=F)

t2sdtL.4 <- cbind(t2sdtL.3, indices)
# t2sdtL.4 %>% kable(digits=2)

t2sdtL.G <- t2sdtL.4 %>%
  dplyr::select(SN, cCong, mPerf, Hit, FA, dprime, aprime, beta, c) %>%
  group_by(cCong) %>%
  summarise(m.p = mean(mPerf), m.d = mean(dprime), m.a = mean(aprime), m.b = mean(beta), m.c = mean(c)) %>% ungroup()
t2sdtL.G %>% kable(digits=2)

cCong	m.p	m.d	m.a	m.b	m.c
v	0.51	1.62	0.85	2.34	0.50
nv	0.47	1.51	0.84	2.37	0.55

# hit - fa
t2sdtW.1 <- t2sdtL.4 %>% dplyr::select(SN, cCong, mPerf) %>% spread(key= cCong, value=mPerf)
# t2sdtW.1 %>% kable(digits=2)

t2sdtG.1 <- Rmisc::summarySEwithin(data = t2sdtL.4, measurevar = "mPerf", idvar = "SN", withinvars = c("cCong"))
# t2sdtG.1 %>% kable(digits=2)

t2.plot2 <- ggplot(data=t2sdtL.4, aes(x=cCong, y=mPerf, fill=cCong)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = 1,
               width = 0.8, show.legend = FALSE) + # , colour="white", size = 1, show.legend = FALSE, colour="black"
  geom_pointrange(data=t2sdtG.1, aes(x = cCong, y=mPerf, ymin = mPerf-se, ymax = mPerf+se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) + #"darkblue"
  geom_segment(data=t2sdtW.1, inherit.aes = FALSE,
               aes(x=1, y=v,
                   xend=2, yend=nv),
               color="gray60") +
  geom_point(data=t2sdtL.4, position=position_dodge(0.80), color="gray60", size=2, show.legend = FALSE) +
  scale_x_discrete(labels=c("Violation", "Non-violation")) +
  scale_fill_manual(values =  c("#feb24c", "#91bfdb")) +
  coord_cartesian(ylim = c(0, 1), clip = "on") +
  # scale_y_continuous(breaks = c(0,0.25,0.50,75,100)) +
  labs(x = "Violation Condition", y = "Hit - FA") +
  # ggtitle("Association Memory Accuracy") +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"),
        # plot.title = element_text(hjust = 0.5),
        legend.position=c(0.8, 0.85))

# d prime
t2sdtW.2 <- t2sdtL.4 %>% dplyr::select(SN, cCong, dprime) %>% spread(key= cCong, value=dprime)
# t2sdtW.2 %>% kable(digits=2)

t2sdtG.2 <- Rmisc::summarySEwithin(data = t2sdtL.4, measurevar = "dprime", idvar = "SN", withinvars = c("cCong"))
# t2sdtG.2 %>% kable(digits=2)

t2.plot3 <- ggplot(data=t2sdtL.4, aes(x=cCong, y=dprime, fill=cCong)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = 1,
               width = 0.8, show.legend = FALSE) + # , colour="white", size = 1, show.legend = FALSE, colour="black"
  geom_pointrange(data=t2sdtG.2, aes(x = cCong, y=dprime, ymin = dprime-se, ymax = dprime+se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) + #"darkblue"
  geom_segment(data=t2sdtW.2, inherit.aes = FALSE,
               aes(x=1, y=v,
                   xend=2, yend=nv),
               color="gray60") +
  geom_point(data=t2sdtL.4, position=position_dodge(0.80), color="gray60", size=2, show.legend = FALSE) +
  scale_x_discrete(labels=c("Violation", "Non-violation")) +
  scale_fill_manual(values =  c("#feb24c", "#91bfdb")) +
  coord_cartesian(ylim = c(0, 4), clip = "on") +
  # scale_y_continuous(breaks = c(0,0.25,0.50,75,100)) +
  labs(x = "Violation Condition", y = "d prime") +
  # ggtitle("Association Memory Accuracy") +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"),
        # plot.title = element_text(hjust = 0.5),
        legend.position=c(0.8, 0.85))
# t2.plot3
t2.plot2.all <-ggarrange(t2.plot2, t2.plot3, ncol = 2)
t2.plot2.all

# criterion
t2sdtW.3 <- t2sdtL.4 %>% dplyr::select(SN, cCong, beta) %>% spread(key= cCong, value=beta)
# t2sdtW.3 %>% kable(digits=2)

t2sdtG.3 <- Rmisc::summarySEwithin(data = t2sdtL.4, measurevar = "beta", idvar = "SN", withinvars = c("cCong"))
# t2sdtG.3 %>% kable(digits=2)

t2.plot4 <- ggplot(data=t2sdtL.4, aes(x=cCong, y=beta, fill=cCong)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = 1,
               width = 0.8, show.legend = FALSE) + # , colour="white", size = 1, show.legend = FALSE, colour="black"
  geom_pointrange(data=t2sdtG.3, aes(x = cCong, y=beta, ymin = beta-se, ymax = beta+se),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) + #"darkblue"
  geom_segment(data=t2sdtW.3, inherit.aes = FALSE,
               aes(x=1, y=v,
                   xend=2, yend=nv),
               color="gray60") +
  geom_point(data=t2sdtL.4, position=position_dodge(0.80), color="gray60", size=2, show.legend = FALSE) +
  scale_x_discrete(labels=c("Violation", "Non-violation")) +
  scale_fill_manual(values =  c("#feb24c", "#91bfdb")) +
  coord_cartesian(ylim = c(0, 4), clip = "on") +
  # scale_y_continuous(breaks = c(0,0.25,0.50,75,100)) +
  labs(x = "Violation Condition", y = "Beta (response bias)") +
  # ggtitle("Association Memory Accuracy") +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"),
        # plot.title = element_text(hjust = 0.5),
        legend.position=c(0.8, 0.85))
t2.plot4

t2accL.1 <- t2sdtL.1.1 # %>% filter(SN!=4) 
t2.dpr.ttest <- t.test(Hit ~ cCong, data=t2accL.1,
                      alternative = c("two.sided"),
                      paired=T, var.equal=F,
                      conf.level=.95)
t2.dpr.ttest

## 
##  Paired t-test
## 
## data:  Hit by cCong
## t = 2.7136, df = 31, p-value = 0.01077
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.009315413 0.065684587
## sample estimates:
## mean of the differences 
##                  0.0375

RM-ANOVA 1

t2accL.1 <- t2sdtL.1.1
t2.acc.aov1 <- aov_ez(id="SN", dv = "Hit", data = t2accL.1, within = c("cCong"))
# test_sphericity(t2.acc.aov1)
# summary(t2.acc.aov1)
# anova(t2.acc.aov1, es = "pes")
nice(t2.acc.aov1, es="pes") %>% kable(digits=3)

Effect	df	MSE	F	pes	p.value
cCong	1, 31	0.00	7.36 *	.192	.011

RM ANOVA 1 - post-hoc

t2.acc.m1 <- emmeans(t2.acc.aov1, pairwise ~ cCong) # adjust="bon" , "tukey"
t2.acc.m1$contrasts %>% kable(digits=3)

contrast	estimate	SE	df	t.ratio	p.value
v - nv	0.038	0.014	31	2.714	0.011

plot(t2.acc.m1, horizontal = FALSE, comparisons = T)

2.2.3 RT Analysis

tt2 <- t2 %>% filter(Corr == 1, cType=="old") 
tt3 <- tt2 %>% filter(RT > 200 & RT < 10000) %>%
  group_by(SN) %>% # grouping by participants
  nest() %>%
  mutate(lbound = map(data, ~mean(.$RT)-3*sd(.$RT)),
         ubound = map(data, ~mean(.$RT)+3*sd(.$RT))) %>% # make new data (3sd cut)
  unnest(c(lbound, ubound))%>%
  unnest(data) %>%
  mutate(Outlier = (RT < lbound)|(RT > ubound)) %>% # set outlier
  filter(Outlier == FALSE) %>% # filtering outlier
  ungroup()

100-100*(nrow(tt3)/nrow(tt2))
## [1] 2.111324

# subject-level, long format (SN/Btw/Block)
t2rtL <- tt2 %>% group_by(SN, cCong) %>%
  dplyr::summarise(RT = mean(RT)) %>%
  ungroup()
## `summarise()` has grouped output by 'SN'. You can override using the `.groups` argument.
# t2rtL %>% kable(digits=2)

# subject-level, wide format (SN/Btw/Block)
t2rtW1 <- t2rtL %>% spread(key=cCong, value = RT)
# t2rtW1 %>% kable(digits=2)

# summary table: grand mean (eyerep/locrep)
t2rtG <- t2rtL %>% group_by(cCong) %>%
  summarise(RT.m = mean(RT), RT.sd = sd(RT)) %>%
  ungroup()
t2rtG$RT.se <- Rmisc::summarySEwithin(data = t2rtL, measurevar = "RT", 
                                        idvar = "SN", withinvars = c("cCong"))$se
t2rtG$RT.ci <- Rmisc::summarySEwithin(data = t2rtL, measurevar = "RT", 
                                        idvar = "SN", withinvars = c("cCong"))$ci
t2rtG <- t2rtG %>% 
  mutate(lower.ci = RT.m-RT.ci,
         upper.ci = RT.m+RT.ci)
# for between-subject design. check help(summarySE)
t2rtG %>% kable(digits=2)

cCong	RT.m	RT.sd	RT.se	RT.ci	lower.ci	upper.ci
v	1573.0	397.24	27.93	56.97	1516.03	1629.98
nv	1591.6	403.22	27.93	56.97	1534.63	1648.58

t2.rt.plot1 <- ggplot(data=t2rtL, aes(x=cCong, y=RT, fill=cCong)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, 
               width = 0.8,size = 1.02, show.legend = FALSE ) + # , show.legend = FALSE  colour="black", 
  geom_point(data=t2rtL, aes(x=cCong, y=RT, fill=cCong), 
             position=position_dodge(width = 0.8), 
             size=2, show.legend = FALSE, color="gray60") + # , 
  geom_segment(data=t2rtW1, inherit.aes = FALSE,
               aes(x=1, y=v,
                   xend=2, yend=nv),
               color="gray60") +
  geom_pointrange(data=t2rtG, aes(x = cCong, y=RT.m, ymin = lower.ci, ymax=upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("Violation", "Non-violation")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4")) +# c("#6798D4", "#D87B49", "#4FA77C"), 
  # scale_fill_manual(values = c("#feb24c", "#91bfdb"), #6798D4 #D87B49 #4FA77C
  # labels = c("Violation", "Non-violation")) +
  coord_cartesian(ylim = c(0, 3000), clip = "on") +
  # scale_y_continuous(breaks = c(0,25,50,75,100)) +
  labs(x = "Stimuli Type", y = "Response Time (ms)", fill ="Violation Condition") +
  # ggtitle("Association Memory Accuracy") +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"), 
        # plot.title = element_text(hjust = 0.5),
        legend.position="top")
# legend.position=c(0.8, 0.90))
t2.rt.plot1

Paired T-Test

t2.rt.ttest <- t.test(RT ~ cCong, data=t2rtL,
                       alternative = c("two.sided"),
                       paired=T, var.equal=F,
                       conf.level=.95)
t2.rt.ttest

## 
##  Paired t-test
## 
## data:  RT by cCong
## t = -0.40773, df = 31, p-value = 0.6863
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -111.63683   74.43789
## sample estimates:
## mean of the differences 
##               -18.59947

---

2.3 Relating priming effect to memory performance

# categorization
t1 <- read.csv("data/reactMEM_p1_all.csv", header = T)
t1$SN = factor(t1$SN)
t1$Block = factor(t1$Block)
t1$cCong = factor(t1$cCong, levels=c(1,2), labels=c("v","nv")) 
t1$cnRep = factor(t1$cnRep, levels=c(1,2,3,5,4), labels=c(1,2,3,4,0))
t1$cnPair = factor(t1$cnPair)
t1$cIdt = factor(t1$cIdt)
t1$cCat = factor(t1$cCat)
t1$Im_Name = factor(t1$Im_Name)
t1$Im_Idx = factor(t1$Im_Idx)
t1$RT <- t1$RT*1000;
t1$Corr <- as.numeric(t1$Corr==1)

# recognition
t2 <- read.csv("data/reactMEM_p2_all.csv", header = T)
t2$SN = factor(t2$SN)
t2$cType = factor(t2$cType, levels=c(1,2), labels=c("old","new"))
t2$cCong = factor(t2$cCong, levels=c(1,2,999), labels=c("v","nv","new")) 
t2$cIdt = factor(t2$cIdt)
t2$dIm_Name = factor(t2$dIm_Name)
t2$Resp <- factor(t2$Resp);
t2$mResp <- factor(t2$Resp, levels=c(1,2,3,4), labels=c("new", "new", "old", "old"))
t2$RT <- t2$RT*1000;
t2$Corr <- as.numeric(t2$Corr==1)
t2$Conf <- as.double(t2$Conf)
t2$p_resp = factor(t2$cType, levels=c('old','new'), labels=c("p.old","p.new"))
t2$r_resp = factor(t2$mResp, levels=c('old','new'), labels=c("r.old", "r.new"))

# RT preprocessing
t3.1 <- t1 %>% filter(Corr == 1) %>% # filter(RT > 200 & RT < 10000) %>%
  group_by(SN) %>% # grouping by participants
  nest() %>%
  mutate(lbound = map(data, ~mean(.$RT)-3*sd(.$RT)),
         ubound = map(data, ~mean(.$RT)+3*sd(.$RT))) %>% # make new data (3sd cut)
  unnest(c(lbound, ubound))%>%
  unnest(data) %>%
  mutate(Outlier = (RT < lbound)|(RT > ubound)) %>% # set outlier
  filter(Outlier == FALSE) %>% # filtering outlier
  ungroup()
t3.1 <- t3.1 %>% filter(cnRep != 1, cnRep != 2)
t3.2 <- t3.1 %>% filter(cnPair == 2) %>% 
  select(SN, cCong, cnRep, cIdt, Im_Idx, RT)
t3.3 <- t3.2 %>% group_by(SN) %>% 
  spread(key = cnRep, value = RT) %>% ungroup()
t3.4 <- t3.3 %>% group_by(SN, cIdt) %>% mutate(prime = `4` - `3`)

# Acc preprocessing
t4.1 <- t2 %>% filter(cType == 'old') %>% 
  select (SN, cCong, cIdt, dIm_Name, Corr, RT, Conf)

t5 <- merge(t4.1, t3.4, key=c(SN,cCong,cIdt))

2.3.1 Logistic Regression

p1.glm <- t5 #  %>% filter(SN!=4)
p1.glm <- na.omit(p1.glm)
table(p1.glm$SN,p1.glm$cCong)

##     
##       v nv new
##   1  46 57   0
##   2  53 53   0
##   3  55 54   0
##   4  50 53   0
##   5  56 59   0
##   6  50 56   0
##   7  45 51   0
##   8  56 54   0
##   9  53 55   0
##   10 46 46   0
##   11 49 53   0
##   12 33 39   0
##   13 54 54   0
##   14 46 47   0
##   15 55 53   0
##   16 51 53   0
##   17 53 53   0
##   18 46 52   0
##   19 55 57   0
##   20 58 54   0
##   21 57 58   0
##   22 55 55   0
##   23 53 52   0
##   24 48 52   0
##   25 54 58   0
##   26 43 48   0
##   27 56 55   0
##   28 57 52   0
##   29 56 52   0
##   30 52 52   0
##   31 50 54   0
##   32 58 56   0

table(p1.glm$Corr,p1.glm$cCong)

##    
##       v  nv new
##   0 724 811   0
##   1 925 886   0

glimpse(p1.glm, width=70)

## Rows: 3,346
## Columns: 12
## $ SN       <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ cCong    <fct> nv, nv, nv, nv, nv, nv, nv, nv, nv, nv, nv, nv, nv,…
## $ cIdt     <fct> 10, 100, 101, 106, 107, 108, 114, 115, 119, 120, 17…
## $ dIm_Name <fct> 144_2_out.jpg, 136_2_out.jpg, 155_2_out.jpg, 101_2_…
## $ Corr     <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, …
## $ RT       <dbl> 1211.4, 969.5, 982.7, 2727.3, 1109.6, 1099.3, 1000.…
## $ Conf     <dbl> 1, 1, 2, 1, 2, 2, 2, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, …
## $ Im_Idx   <fct> 144, 136, 155, 101, 194, 132, 3, 172, 32, 14, 22, 1…
## $ `3`      <dbl> 529.5, 568.3, 538.0, 479.2, 536.1, 506.8, 478.5, 51…
## $ `4`      <dbl> 479.1, 531.8, 670.9, 494.6, 526.4, 467.2, 457.5, 43…
## $ `0`      <dbl> 481.6, 538.8, 470.1, 521.0, 538.9, 536.4, 587.7, 44…
## $ prime    <dbl> -50.4, -36.5, 132.9, 15.4, -9.7, -39.6, -21.0, -80.…

p1.glm.v <- p1.glm %>% filter(cCong == "v")
p1.glm1 <- glm(Corr ~ prime, data=p1.glm.v, family=binomial()) 
# summary(p1.glm1)  %>% kable(digits=2)
# sjPlot::tab_model(p1.glm1)

p1.glm1.2 <- glm(Corr ~ `0`, data=p1.glm.v, family=binomial())
# summary(p1.glm1.2) # display results
# anova(p1.glm1.2)
# sjPlot::tab_model(p1.glm1.2)

p1.glm.v$prime2 <- p1.glm.v$prime^2
p1.glm1.1 <- glm(Corr ~ prime + prime2, data=p1.glm.v, family=binomial()) 
# summary(p1.glm1.1) # display results
# anova(p1.glm1.1)
# sjPlot::tab_model(p1.glm1.1)
sjPlot::tab_model(p1.glm1)

	Corr
Predictors	Odds Ratios	CI	p
(Intercept)	1.27	1.16 – 1.41	<0.001
prime	1.00	1.00 – 1.00	0.240
Observations	1649
R2 Tjur	0.001

sjPlot::tab_model(p1.glm1.1)

	Corr
Predictors	Odds Ratios	CI	p
(Intercept)	1.30	1.17 – 1.44	<0.001
prime	1.00	1.00 – 1.00	0.218
prime2	1.00	1.00 – 1.00	0.258
Observations	1649
R2 Tjur	0.002

sjPlot::tab_model(p1.glm1.2)

	Corr
Predictors	Odds Ratios	CI	p
(Intercept)	1.63	1.20 – 2.22	0.002
0	1.00	1.00 – 1.00	0.099
Observations	1649
R2 Tjur	0.002

p3_logreg.v <- ggplot(p1.glm.v, aes(x=prime, y=Corr)) + 
  geom_point(alpha=.5, size=2) +
  # ggbeeswarm::geom_quasirandom(dodge.width = 1, color = "blue", size = 1, alpha = 0.2,
                               # show.legend = FALSE) +
  stat_smooth(method="glm", se=T, method.args = list(family=binomial)) + 
  labs(x = "Diff Score (final B rt - 3rd B rt)", y = "Correctness") +
  scale_y_continuous(breaks = c(0, 0.5, 1)) +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"), 
        # plot.title = element_text(hjust = 0.5),
        legend.position=c(0.8, 0.85))
# p3_logreg.v

p3_logreg.v.2 <- ggplot(p1.glm.v, aes(x=`0`, y=Corr)) + 
  geom_point(alpha=.5, size=2) +
  # ggbeeswarm::geom_quasirandom(dodge.width = 1, color = "blue", size = 1, alpha = 0.2,
  # show.legend = FALSE) +
  stat_smooth(method="glm", se=T, method.args = list(family=binomial)) + 
  labs(x = "B` RT", y = "Correctness") +
  scale_y_continuous(breaks = c(0, 0.5, 1)) +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"), 
        # plot.title = element_text(hjust = 0.5),
        legend.position=c(0.8, 0.85))


p3_logreg.plot <- ggarrange(p3_logreg.v, p3_logreg.v.2, ncol = 2,
          labels=c("A) Violation, Prime", "B) Violation, B`"))

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

p3_logreg.plot

p1.glm.nv <- p1.glm %>% filter(cCong == "nv")
p1.glm.nv <- p1.glm.nv %>% group_by(SN) %>% mutate(prime_z = (prime-mean(prime))/sd(prime)) %>% ungroup()

p1.glm2 <- glm(Corr ~ prime, data=p1.glm.nv, family=binomial()) 
# summary(p1.glm2) # display results
# anova(p1.glm2)

p1.glm.nv$prime2 <- p1.glm.nv$prime^2
p1.glm2.1 <- glm(Corr ~ prime + prime2, data=p1.glm.nv, family=binomial()) 
# summary(p1.glm2.1) # display results
# anova(p1.glm2.1)

p1.glm2.2 <- glm(Corr ~ `0`, data=p1.glm.nv, family=binomial())
# summary(p1.glm2.2) # display results
# anova(p1.glm2.2)

sjPlot::tab_model(p1.glm2)

	Corr
Predictors	Odds Ratios	CI	p
(Intercept)	1.09	0.99 – 1.20	0.067
prime	1.00	1.00 – 1.00	0.753
Observations	1697
R2 Tjur	0.000

sjPlot::tab_model(p1.glm2.1)

	Corr
Predictors	Odds Ratios	CI	p
(Intercept)	1.09	0.98 – 1.21	0.096
prime	1.00	1.00 – 1.00	0.764
prime2	1.00	1.00 – 1.00	0.872
Observations	1697
R2 Tjur	0.000

sjPlot::tab_model(p1.glm2.2)

	Corr
Predictors	Odds Ratios	CI	p
(Intercept)	1.61	1.20 – 2.17	0.002
0	1.00	1.00 – 1.00	0.007
Observations	1697
R2 Tjur	0.004

# hist(residuals(p1.glm1, type="deviance"))
# qqnorm(residuals(p1.glm1, type="deviance"))

# plotting
newdat.nv <- data.frame(prime=seq(min(p1.glm.nv$prime), max(p1.glm.nv$prime),len=dim(p1.glm.nv)[1]))
newdat.nv$Corr = predict(p1.glm2, newdata=p1.glm.nv, type="response")

p3_logreg.nv <- ggplot(p1.glm.nv, aes(x=prime, y=Corr)) + 
  geom_point(alpha=.5, size=2) +
  # ggbeeswarm::geom_quasirandom(dodge.width = 1, color = "blue", size = 1, alpha = 0.2,
  # show.legend = FALSE) +
  stat_smooth(method=glm, se=T, method.args = list(family=binomial)) + 
  labs(x = "Diff Score (final B rt - 3rd B rt)", y = "Correctness") +
  scale_y_continuous(breaks = c(0, 0.5, 1)) +
  # ggtitle("Non-violation Condition") +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"), 
        # plot.title = element_text(hjust = 0.5),
        legend.position=c(0.8, 0.85))
# p3_logreg.nv

p3_logreg.nv.2 <- ggplot(p1.glm.nv, aes(x=`0`, y=Corr)) + 
  geom_point(alpha=.5, size=2) +
  # ggbeeswarm::geom_quasirandom(dodge.width = 1, color = "blue", size = 1, alpha = 0.2,
  # show.legend = FALSE) +
  stat_smooth(method="glm", se=T, method.args = list(family=binomial)) + 
  labs(x = "B` RT", y = "Correctness") +
  scale_y_continuous(breaks = c(0, 0.5, 1)) +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"), 
        # plot.title = element_text(hjust = 0.5),
        legend.position=c(0.8, 0.85))
# p3_logreg.nv.2

p3_logreg.plot <- ggarrange(p3_logreg.nv, p3_logreg.nv.2, ncol = 2,
          labels=c("A) non-Violation, Prime", "B) non-Violation, B`"))

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

p3_logreg.plot

2.3.2 Generalized Linear Mixed Effect Model

p1.glm.v <- p1.glm.v %>% group_by(SN) %>% 
  mutate(prime_z = (prime-mean(prime))/sd(prime)) %>% ungroup()

(nc <- detectCores())

## [1] 8

cl <- makeCluster(rep("localhost", nc))
p1.glm <- p1.glm %>% group_by(SN) %>%
  mutate(prime_z_a = (prime-mean(prime))/sd(prime)) %>% ungroup()
p1.glm.mixed <- afex::mixed(Corr ~  cCong*prime + (1|SN) + (1|Im_Idx),
                              data = p1.glm,
                              family = binomial(link="logit"),
                              method = "LRT", cl = cl,
                              expand_re = T,
                              control = lmerControl(optimizer = "bobyqa",
                                                    optCtrl = list(maxfun = 1e7)))

## Fitting 4 (g)lmer() models.

p1.glm.mixed.v <- afex::mixed(Corr ~ prime_z + (1|SN) + (1|Im_Idx),
                              data = p1.glm.v,
                              family = binomial(link="logit"),
                              method = "LRT", cl = cl,
                              expand_re = T,
                              control = lmerControl(optimizer = "bobyqa",
                                                    optCtrl = list(maxfun = 1e7)))

## Fitting 2 (g)lmer() models.

p1.glm.mixed.nv <- afex::mixed(Corr ~ prime_z + (1|SN) + (1|Im_Idx),
                               data = p1.glm.nv,
                               family = binomial(link="logit"),
                               method = "LRT", cl = cl,
                               expand_re = T,
                               control = lmerControl(optimizer = "bobyqa",
                                                     optCtrl = list(maxfun = 1e7)))

## Fitting 2 (g)lmer() models.

stopCluster(cl)

sjPlot::tab_model(p1.glm.mixed)

	Corr
Predictors	Estimates	CI	p
(Intercept)	0.28	0.00 – 0.56	0.046
cCong [nv]	-0.17	-0.32 – -0.02	0.031
prime	0.00	-0.00 – 0.00	0.435
cCong [nv] * prime	-0.00	-0.00 – 0.00	0.740
Random Effects
σ2	3.29
τ00 Im_Idx	0.30
τ00 SN	0.50
ICC	0.20
N SN	32
N Im_Idx	240
Marginal R2 / Conditional R2	0.002 / 0.197

sjPlot::tab_model(p1.glm.mixed.v)

	Corr
Predictors	Estimates	CI	p
(Intercept)	0.31	0.04 – 0.58	0.027
prime_z	0.05	-0.06 – 0.16	0.387
Random Effects
σ2	3.29
τ00 Im_Idx	0.36
τ00 SN	0.47
ICC	0.20
N SN	32
N Im_Idx	240
Marginal R2 / Conditional R2	0.001 / 0.202

sjPlot::tab_model(p1.glm.mixed.nv)

	Corr
Predictors	Estimates	CI	p
(Intercept)	0.11	-0.15 – 0.38	0.409
prime_z	0.06	-0.05 – 0.16	0.271
Random Effects
σ2	3.29
τ00 Im_Idx	0.17
τ00 SN	0.48
ICC	0.17
N SN	32
N Im_Idx	240
Marginal R2 / Conditional R2	0.001 / 0.166

2.3.3 Binning and ANOVA

p1.bin <- p1.glm  %>%  group_by(SN, cCong) %>% # grouping by participants
  mutate(med = median(prime)) %>% # make new data (3sd cut)
  mutate(prime_bin = (prime >= med)) %>% # set outlier
  ungroup()

p1.bin$prime_bin <- as.numeric(p1.bin$prime_bin == T)
p1.bin$prime_bin <- p1.bin$prime_bin+1
p1.bin$prime_bin <- factor(p1.bin$prime_bin, levels=c(1,2), labels = c("low", "high"))
glimpse(p1.bin)

## Rows: 3,346
## Columns: 15
## $ SN        <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ cCong     <fct> nv, nv, nv, nv, nv, nv, nv, nv, nv, nv, nv, nv, nv, nv, nv, …
## $ cIdt      <fct> 10, 100, 101, 106, 107, 108, 114, 115, 119, 120, 17, 18, 19,…
## $ dIm_Name  <fct> 144_2_out.jpg, 136_2_out.jpg, 155_2_out.jpg, 101_2_in.jpg, 1…
## $ Corr      <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, …
## $ RT        <dbl> 1211.4, 969.5, 982.7, 2727.3, 1109.6, 1099.3, 1000.0, 799.1,…
## $ Conf      <dbl> 1, 1, 2, 1, 2, 2, 2, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2, …
## $ Im_Idx    <fct> 144, 136, 155, 101, 194, 132, 3, 172, 32, 14, 22, 141, 5, 15…
## $ `3`       <dbl> 529.5, 568.3, 538.0, 479.2, 536.1, 506.8, 478.5, 514.4, 538.…
## $ `4`       <dbl> 479.1, 531.8, 670.9, 494.6, 526.4, 467.2, 457.5, 433.8, 509.…
## $ `0`       <dbl> 481.6, 538.8, 470.1, 521.0, 538.9, 536.4, 587.7, 448.5, 482.…
## $ prime     <dbl> -50.4, -36.5, 132.9, 15.4, -9.7, -39.6, -21.0, -80.6, -28.9,…
## $ prime_z_a <dbl> -0.55233883, -0.40079548, 1.44607094, 0.16503904, -0.1086111…
## $ med       <dbl> 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, …
## $ prime_bin <fct> low, low, high, high, low, low, low, low, low, high, low, lo…

# subject-level, long format (SN/Btw/Block)
t5rtL <- p1.bin %>% group_by(SN, cCong, prime_bin) %>%
  dplyr::summarise(p2Acc = mean(Corr)*100) %>%
  ungroup()
# t5rtL %>% kable(digits=2)

# subject-level, wide format (SN/Btw/Block)
t5rtW1 <- t5rtL %>% spread(key=prime_bin, value = p2Acc)
# t5rtW1 %>% kable(digits=2)

# summary table: grand mean (eyerep/locrep)
t5rtG <- t5rtL %>% group_by(cCong, prime_bin) %>%
  summarise(p2Acc.m = mean(p2Acc), p2Acc.sd = sd(p2Acc)) %>%
  ungroup()
t5rtG$p2Acc.se <- Rmisc::summarySEwithin(data = t5rtL, measurevar = "p2Acc", 
                                      idvar = "SN", withinvars = c("cCong","prime_bin"))$se
t5rtG$p2Acc.ci <- Rmisc::summarySEwithin(data = t5rtL, measurevar = "p2Acc", 
                                      idvar = "SN", withinvars = c("cCong","prime_bin"))$ci
t5rtG <- t5rtG %>% 
  mutate(lower.ci = p2Acc.m-p2Acc.ci,
         upper.ci = p2Acc.m+p2Acc.ci)
# for between-subject design. check help(summarySE)
t5rtG %>% kable(digits=2)

cCong	prime_bin	p2Acc.m	p2Acc.sd	p2Acc.se	p2Acc.ci	lower.ci	upper.ci
v	low	55.71	18.44	1.78	3.64	52.07	59.35
v	high	56.78	16.06	1.22	2.50	54.28	59.27
nv	low	51.99	18.71	1.56	3.18	48.81	55.17
nv	high	52.59	17.92	1.55	3.15	49.43	55.74

t5.plot1 <- ggplot(data=t5rtL, aes(x=cCong, y=p2Acc, fill=prime_bin)) +
  stat_summary(fun = mean, geom = "bar", position="dodge", na.rm = TRUE, alpha = .9, 
               width = 0.8,size = 1.02) + # , show.legend = FALSE  colour="black", 
  geom_point(data=t5rtL, aes(x=cCong, y=p2Acc, fill=prime_bin), 
             position=position_dodge(width = 0.8), 
             size=2, show.legend = FALSE, color="gray60") + # , 
  geom_segment(data=filter(t5rtW1, cCong == "v"), inherit.aes = FALSE,
               aes(x=0.8, y=filter(t5rtW1, cCong == "v")$low,
                   xend=1.2, yend=filter(t5rtW1, cCong == "v")$high),
               color="gray60") +
  
  geom_segment(data=filter(t5rtW1, cCong == "nv"), inherit.aes = FALSE,
               aes(x=1.8, y=filter(t5rtW1, cCong == "nv")$low,
                   xend=2.2, yend=filter(t5rtW1, cCong == "nv")$high),
               color="gray60") +
  
  geom_pointrange(data=t5rtG, aes(x = cCong, y=p2Acc.m, ymin = lower.ci, ymax = upper.ci),
                  position = position_dodge(0.80), color = "darkred", size = 1, show.legend = FALSE) +
  scale_x_discrete(labels=c("Violation","Non-violation")) +
  scale_fill_manual(values = c("#D2E6BD", "#4AA7B4"),# c("#6798D4", "#D87B49", "#4FA77C"), 
                    labels = c("Low", "High")) +
  # scale_fill_manual(values = c("#feb24c", "#91bfdb"), #6798D4 #D87B49 #4FA77C
  # labels = c("Violation", "Non-violation")) +
  coord_cartesian(ylim = c(0, 120), clip = "on") +
  scale_y_continuous(breaks = c(0,25,50,75,100)) +
  labs(x = "Violation Condition", y = "p2 Recognition Hit (%)", fill ="p1 4th-3rd B rt") +
  # ggtitle("Association Memory Accuracy") +
  theme_bw(base_size = 15) +
  theme(axis.title = element_text(face = "bold", size = 16, color = "black"),
        axis.text = element_text(face = "plain", hjust = 0.5, size = 15, color = "black"),
        axis.line=element_line(),
        strip.text.x = element_text(face = "plain", size = 15, color = "black"),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing=unit(1, "lines"),
        # aspect.ratio = 0.2,
        plot.margin = margin(0.3, 0.3, 0.3, 0.3, "cm"), 
        # plot.title = element_text(hjust = 0.5),
        legend.position="top")
# legend.position=c(0.8, 0.90))
t5.plot1

RM-ANOVA

t5.acc.aov1 <- aov_ez(id="SN", dv = "p2Acc", data = t5rtL, within = c("cCong", "prime_bin"))
# test_sphericity(t2.acc.aov1)
# summary(t5.acc.aov1)
# anova(t5.acc.aov1, es = "pes")
nice(t5.acc.aov1, es="pes") %>% kable(digits=3)

Effect	df	MSE	F	pes	p.value
cCong	1, 31	76.33	6.55 *	.174	.016
prime_bin	1, 31	84.81	0.26	.008	.614
cCong:prime_bin	1, 31	92.33	0.02	<.001	.890

RM ANOVA - post-hoc

t5.acc.m1 <- emmeans(t5.acc.aov1, pairwise ~ prime_bin | cCong) # adjust="bon" , "tukey"
t5.acc.m1$contrasts %>% kable(digits=3)

contrast	cCong	estimate	SE	df	t.ratio	p.value
low - high	v	-1.065	2.328	31	-0.457	0.651
low - high	nv	-0.593	2.378	31	-0.249	0.805

plot(t5.acc.m1, horizontal = FALSE, comparisons = T)

Paired T-Test

Violation

t5rtL.v <- t5rtL %>% filter(cCong=="v")
t5.acc.ttest <- t.test(p2Acc ~ prime_bin, data=t5rtL.v,
                       alternative = c("two.sided"),
                       paired=T, var.equal=F,
                       conf.level=.95)
t5.acc.ttest

## 
##  Paired t-test
## 
## data:  p2Acc by prime_bin
## t = -0.45735, df = 31, p-value = 0.6506
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.812061  3.682881
## sample estimates:
## mean of the differences 
##                -1.06459

Non-Violation

t5rtL.nv <- t5rtL %>% filter(cCong=="nv")

t5.acc.ttest <- t.test(p2Acc ~ prime_bin, data=t5rtL.nv,
                       alternative = c("two.sided"),
                       paired=T, var.equal=F,
                       conf.level=.95)
t5.acc.ttest

## 
##  Paired t-test
## 
## data:  p2Acc by prime_bin
## t = -0.24923, df = 31, p-value = 0.8048
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.441678  4.256527
## sample estimates:
## mean of the differences 
##              -0.5925753